Code
pacman::p_load(sf, tidyverse, tmap, dplyr,
spatstat, spdep,
lubridate, leaflet,
plotly, DT, viridis,
ggplot2, sfdep, janitor, shiny, rnaturalearth)Andrea Yeo
March 22, 2025
March 27, 2025
For this take-home exercise, each team member have to choose one module from our proposed Shiny application and complete the following tasks:
Our project utilizes open-source data from the World Happiness Report – Data Sharing covering the years 2015 to 2025. Detailed information on the data cleaning and preparation process is available on our group’s Netlify site, under the “Data Preparation” tab: Group Visual Analytic Application.
Specifically, our project focuses on the visualisation of global happiness trends. I will be exploring the modules related to geospatial and aspatial analysis to explore both spatial patterns and non-spatial factors influencing happiness scores across countries. For the prototyping of these modules, I will be using the cleaned_data dataset prepared after the data processing phase.
The below R packages will be used in this exercise and for the Shiny application
glimpse(): provides a transposed overview of a dataset, showing variables and their types in a concise format.head(): displays the first few rows of a dataset (default is 6 rows) to give a quick preview of the data.summary(): generates a statistical summary of each variable, including measures like mean, median, and range for numeric data.duplicated():returns a logical vector indicating which elements or rows in a vector or data frame are duplicates.colSums(is.na()): counts the number of missing values (NA) in each column of the data frame.str(): use str() to display the column names, data types, and a preview of the data.Rows: 1,656
Columns: 12
$ year <int> 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2…
$ country <chr> "Afghanistan", "Afghanistan", "Afghanistan", "Afg…
$ ladder_score <dbl> 3.5750, 3.3600, 3.7940, 3.6320, 3.2030, 2.5669, 2…
$ economy_score <dbl> 0.31982, 0.38227, 0.40148, 0.33200, 0.35000, 0.30…
$ social_score <dbl> 0.30285, 0.11037, 0.58154, 0.53700, 0.51700, 0.35…
$ lifeexpectancy_score <dbl> 0.30335, 0.17344, 0.18075, 0.25500, 0.36100, 0.26…
$ freedom_score <dbl> 0.23414, 0.16430, 0.10618, 0.08500, 0.00000, 0.00…
$ generosity_score <dbl> 0.36510, 0.31268, 0.31187, 0.19100, 0.15800, 0.13…
$ corrperception_score <dbl> 0.09719, 0.07112, 0.06116, 0.03600, 0.02500, 0.00…
$ residual_score <dbl> 1.95210, 2.14558, 2.15080, 2.19600, 1.79200, 1.50…
$ rank <int> 153, 154, 141, 145, 154, 153, 149, 146, 137, 143,…
$ region <chr> "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "…
year country ladder_score economy_score social_score lifeexpectancy_score
1 2014 Afghanistan 3.5750 0.31982 0.30285 0.30335
2 2015 Afghanistan 3.3600 0.38227 0.11037 0.17344
3 2016 Afghanistan 3.7940 0.40148 0.58154 0.18075
4 2017 Afghanistan 3.6320 0.33200 0.53700 0.25500
5 2018 Afghanistan 3.2030 0.35000 0.51700 0.36100
6 2019 Afghanistan 2.5669 0.30071 0.35643 0.26605
freedom_score generosity_score corrperception_score residual_score rank
1 0.23414 0.36510 0.09719 1.95210 153
2 0.16430 0.31268 0.07112 2.14558 154
3 0.10618 0.31187 0.06116 2.15080 141
4 0.08500 0.19100 0.03600 2.19600 145
5 0.00000 0.15800 0.02500 1.79200 154
6 0.00000 0.13523 0.00123 1.50724 153
region
1 Asia
2 Asia
3 Asia
4 Asia
5 Asia
6 Asia
year country ladder_score economy_score
Min. :2014 Length:1656 Min. :1.364 Min. :0.0000
1st Qu.:2016 Class :character 1st Qu.:4.606 1st Qu.:0.7607
Median :2019 Mode :character Median :5.497 Median :1.0983
Mean :2019 Mean :5.460 Mean :1.0763
3rd Qu.:2022 3rd Qu.:6.303 3rd Qu.:1.3946
Max. :2024 Max. :7.842 Max. :2.2090
social_score lifeexpectancy_score freedom_score generosity_score
Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
1st Qu.:0.8642 1st Qu.:0.4028 1st Qu.:0.3794 1st Qu.:0.1069
Median :1.1122 Median :0.5963 Median :0.4972 Median :0.1700
Mean :1.0782 Mean :0.5756 Mean :0.4914 Mean :0.1846
3rd Qu.:1.3496 3rd Qu.:0.7599 3rd Qu.:0.6110 3rd Qu.:0.2421
Max. :1.8400 Max. :1.1410 Max. :1.0180 Max. :0.8381
corrperception_score residual_score rank region
Min. :0.0000 Min. :-0.110 Min. : 1.00 Length:1656
1st Qu.:0.0590 1st Qu.: 1.553 1st Qu.: 38.00 Class :character
Median :0.1010 Median : 1.922 Median : 76.00 Mode :character
Mean :0.1356 Mean : 1.919 Mean : 75.97
3rd Qu.:0.1726 3rd Qu.: 2.309 3rd Qu.:113.25
Max. :0.5870 Max. : 3.838 Max. :158.00
[1] year country ladder_score
[4] economy_score social_score lifeexpectancy_score
[7] freedom_score generosity_score corrperception_score
[10] residual_score rank region
<0 rows> (or 0-length row.names)
year country ladder_score
0 0 0
economy_score social_score lifeexpectancy_score
0 0 0
freedom_score generosity_score corrperception_score
0 0 0
residual_score rank region
0 0 0
drop_na() function to drop rows where any specified column contains a missing value.'data.frame': 1656 obs. of 12 variables:
$ year : int 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 ...
$ country : chr "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
$ ladder_score : num 3.58 3.36 3.79 3.63 3.2 ...
$ economy_score : num 0.32 0.382 0.401 0.332 0.35 ...
$ social_score : num 0.303 0.11 0.582 0.537 0.517 ...
$ lifeexpectancy_score: num 0.303 0.173 0.181 0.255 0.361 ...
$ freedom_score : num 0.234 0.164 0.106 0.085 0 ...
$ generosity_score : num 0.365 0.313 0.312 0.191 0.158 ...
$ corrperception_score: num 0.0972 0.0711 0.0612 0.036 0.025 ...
$ residual_score : num 1.95 2.15 2.15 2.2 1.79 ...
$ rank : int 153 154 141 145 154 153 149 146 137 143 ...
$ region : chr "Asia" "Asia" "Asia" "Asia" ...
The happiness tibble contains 12 attributes, as shown above.
The following preprocessing checks were conducted as part of data preparation:
happiness dataset using glimpse() and str()duplicated() in the datasetcolSums(is.na())To support the geospatial visualisation of global happiness trends, the rnaturalearth package was used to import spatial boundary data of countries. Specifically, the ne_countries() function was called with a medium scale and returned as a simple features (sf) object, stored in the world variable. This spatial dataset provides the geographic outlines of all countries and is essential for merging with the World Happiness Report data to enable choropleth mapping and spatial analysis at the country level.
To prepare the dataset for mapping, the World Happiness data was filtered to include only records from the year 2024. This filtered data (happiness_latest) was then joined with the spatial dataset (world) using the country name as the common key. The resulting dataset, world_happy, combines both geographic boundaries and happiness scores, enabling geospatial visualisation of happiness levels across countries.
The code below plots a choropleth map to visualise the distribution of happiness scores across countries for the year 2024. Using the tmap package in plotting mode, the tm_shape() function defines the spatial data object (world_happy), which contains both geographic boundaries and the associated happiness data. The tm_polygons() function maps the ladder_score variable to a blue color gradient, with darker shades indicating higher happiness scores.
Finally, tm_layout() adds a descriptive title and positions the legend outside the map area to improve clarity and presentation.
The enhanced choropleth map uses the tmap package with improved styling and classification. It applies a quantile-based scheme with five classes and the "YlOrRd" palette to represent ladder_score, where darker shades indicate higher happiness. Additional cartographic elements such as a compass, scale bar, borders, and descriptive credits are included to improve readability and interpretability. The legend is placed on the right, and a classic map style is applied for a cleaner layout.
tmap_style("classic")
tmap_mode("plot")
tm_shape(world_happy) +
tm_polygons(
col = "ladder_score",
palette = "YlOrRd",
style = "quantile",
n = 5,
title = "Happiness Score (2024)"
) +
tm_layout(
main.title = "World Happiness Distribution (2024)",
main.title.position = "center",
main.title.size = 1.2,
legend.outside = TRUE,
legend.position = c("right", "center"),
frame = FALSE
) +
tm_borders(alpha = 0.4, lwd = 0.3) +
tm_compass(type = "8star", size = 2, position = c("left", "top")) +
tm_scale_bar(position = c("left", "bottom")) +
tm_credits("Source: World Happiness Report 2024", position = c("left", "bottom")) +
tm_credits("Note: Darker shades indicate higher happiness; lighter shades indicate lower happiness.",
position = c("left", "bottom"), just = "right", size = 0.6)
| Feature | Basic Version | Enhanced Version |
|---|---|---|
| Classification method | Default(“pretty”) | Custom(“quantile”, 5 bins) |
| Layout | Basic title, default layout | Custom title, centered, larger size |
| Map Borders | None | Added with transparency and thin width |
| Compass/ Scale Bar | Not included | Included |
| Credits/ Source | Not included | Included |
| Frame | Default | Removed |
The code below creates a faceted choropleth map to visualise global happiness scores from 2014 to 2024. Using tmap, each facet represents a different year, allowing for easy comparison of happiness trends across countries. Geospatial data is joined using rnaturalearth.
# Load required libraries
library(tidyverse); library(tmap); library(sf); library(rnaturalearth)
# Set style and plotting mode
tmap_style("white")
tmap_mode("plot")
# Step 1: Filter for all years from 2014 to 2024
years_to_plot <- 2014:2024
happiness_subset <- happiness %>%
filter(year %in% years_to_plot)
# Step 2: Load and join with world spatial data
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy_facet <- world %>%
left_join(happiness_subset, by = c("name" = "country"))
# Step 3: Filter out rows with missing ladder_score
world_happy_facet_clean <- world_happy_facet %>%
filter(!is.na(ladder_score))
# Step 4: Plot small multiple choropleth maps
tm_shape(world_happy_facet_clean) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
style = "quantile",
title = "Happiness Score"
) +
tm_facets(by = "year", ncol = 2) +
tm_layout(
main.title = "Global Happiness Scores (2014–2024)",
main.title.position = "center",
panel.label.size = 2.0,
legend.outside = TRUE,
legend.position = c("center", "top"),
legend.height = 1.5,
legend.width = 1.0,
legend.text.size = 1.5,
legend.title.size = 1.5,
asp = 0,
between.margin = 0,
outer.margins = c(0, 0, 0, 0)
)
The code below produces an animated choropleth map using ggplot2, gganimate, and sf to illustrate the temporal evolution of global happiness scores from 2014 to 2024. Country-level geospatial data from rnaturalearth is merged with World Happiness Report data via a spatial join. transition_time() is used to animate by year, with countries shaded according to their ladder score.
Note: Although the underlying world geometries are static, the animated map may appear to shift or jitter between frames. This occurs because some countries are missing data in certain years, leading to variations in the number of plotted geometries. As a result, gganimate recalculates the plot layout for each frame, even when coord_sf(xlim, ylim, expand = FALSE) is specified. The visual movement is not due to changes in the map itself, but rather to inconsistencies in data availability across years.
library(tidyverse); library(sf); library(gganimate); library(rnaturalearth)
# Load base world map
world <- ne_countries(scale = "medium", returnclass = "sf") %>%
select(name, geometry)
# Load and clean happiness data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year %in% 2014:2024)
# Merge first — keep all geometries
world_happy <- world %>%
left_join(happiness, by = c("name" = "country")) %>%
mutate(year = as.numeric(year))
# Plot
ggplot(world_happy) +
geom_sf(aes(fill = ladder_score)) +
scale_fill_distiller(palette = "YlGnBu", direction = 1, na.value = "lightgray") +
labs(title = "Happiness Score: {frame_time}", fill = "Score") +
transition_time(year) +
ease_aes('linear') +
coord_sf(
xlim = c(-180, 180),
ylim = c(-60, 90),
expand = FALSE
)
This code creates an interactive world map using tmap, showing happiness scores by country for 2024. Each country is colored based on its ladder score, and tooltips display the country name when hovered. The map allows users to zoom and explore happiness data visually.
happiness_latest <- happiness_subset %>%
filter(year == 2024)
world_happy_latest <- world %>%
left_join(happiness_latest, by = c("name" = "country"))
world_happy_latest_clean <- world_happy_latest %>%
filter(!is.na(ladder_score))
tmap_mode("view")
tm_shape(world_happy_latest_clean) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name", # this is what shows in the tooltip
title = "Happiness Score (2024)"
)In this module, I developed an interactive Shiny application to visualize the World Happiness Score by country and year. The key input is a dropdown menu (selectInput) that allows users to choose a specific year from 2014 to 2024. Based on the selected year, the output is a dynamic choropleth map (tmapOutput) showing each country’s happiness score using a color gradient.
Input: A single selectInput() allows the user to choose the Year of interest from 2014 to 2024.
Output: A dynamic choropleth map is rendered using tmap, updating automatically based on the selected year.
UI Components Used: selectInput() for the dropdown menu, sidebarPanel() for input placement, and tmapOutput() in mainPanel() for rendering the interactive map.

# Load libraries
library(shiny); library(tmap); library(tidyverse); library(sf); library(rnaturalearth)
# Prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year))
years_to_plot <- 2014:2024
happiness_subset <- happiness %>%
filter(year %in% years_to_plot)
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
left_join(happiness_subset, by = c("name" = "country"))
# UI
ui <- fluidPage(
titlePanel("🌍 World Happiness Score by Year"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:",
choices = sort(unique(world_happy$year)),
selected = 2024)
),
mainPanel(
tmapOutput("happiness_map")
)
)
)
# Server
server <- function(input, output, session) {
filtered_data <- reactive({
world_happy %>%
filter(year == input$selected_year & !is.na(ladder_score))
})
output$happiness_map <- renderTmap({
tmap_mode("view")
tm_shape(filtered_data()) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name",
popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
title = paste("Happiness Score:", input$selected_year)
)
})
}
# Run the app
shinyApp(ui, server)This code below generates an interactive map using leaflet to visualize 2024 World Happiness Scores. Each country is displayed as a circle located at its geographic center, with the size representing its happiness score. I added tooltips which provide additional insights like economy, life expectancy, freedom, and region, and a legend which helps interpret the happiness scores visually.
library(tidyverse); library(readr); library(leaflet)
library(sf); library(rnaturalearth); library(rnaturalearthdata)
# Load the world happiness data
happiness <- read_csv("data/world_happiness.csv")
# Filter for a specific year (e.g., 2024)
happiness_latest <- happiness %>%
filter(year == 2024)
# Load country polygons and calculate centroids
world_sf <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
select(name, geometry)
# Join happiness data to spatial data
world_happy_sf <- world_sf %>%
left_join(happiness_latest, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
# Get centroids for each country
world_happy_points <- st_centroid(world_happy_sf)
# Extract coordinates for leaflet
world_happy_points_coords <- world_happy_points %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
leaflet(world_happy_points_coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)In this module, I developed an interactive Shiny application that allows users to explore the World Happiness Score by country for the year 2024. The main input is a country selector using selectInput, which enables users to search and select a country from the dropdown list. Once selected, the output is an interactive leaflet map that automatically zooms to the selected country and displays a popup with additional details such as economy, life expectancy, freedom, and region.
Input: A selectInput() enables users to search and select a country from the list of available options for 2024.
Output: A proportional symbol map rendered via leaflet, displaying country-level happiness scores using interactive circle markers with popups.
UI Components Used: selectInput() in the sidebarPanel() for country selection, and leafletOutput() in the mainPanel() for displaying the dynamic map.

library(shiny); library(tidyverse); library(readr); library(leaflet)
library(sf); library(rnaturalearth); library(rnaturalearthdata)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
filter(year == 2024)
world_sf <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
select(name, geometry)
world_happy_sf <- world_sf %>%
left_join(happiness, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
world_happy_points <- st_centroid(world_happy_sf)
world_happy_points_coords <- world_happy_points %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
# UI
ui <- fluidPage(
titlePanel("World Happiness (2024)"),
sidebarLayout(
sidebarPanel(
selectInput("selected_country", "Search Country:",
choices = sort(world_happy_points_coords$name),
selected = NULL)
),
mainPanel(
leafletOutput("map", height = "400px")
)
)
)
# Server
server <- function(input, output, session) {
output$map <- renderLeaflet({
leaflet(world_happy_points_coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
})
observeEvent(input$selected_country, {
selected_data <- world_happy_points_coords %>%
filter(name == input$selected_country)
leafletProxy("map") %>%
setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
clearPopups() %>%
addPopups(
lng = selected_data$lon,
lat = selected_data$lat,
popup = paste0(
"<b>Country:</b> ", selected_data$name, "<br/>",
"<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
"<b>Region:</b> ", selected_data$region
)
)
})
}
# Run the app
shinyApp(ui, server)This module presents an interactive Shiny application with two coordinated maps to explore World Happiness data. The Choropleth Map allows users to select a year and view global happiness patterns using a color gradient. In contrast, the Proportional Symbol Map lets users search for a country, zoom in automatically, and view detailed indicators such as economy, life expectancy, and freedom. Circle size and color reflect the happiness score.
After considering various approaches, I chose to implement both map types as they serve different purposes:
the choropleth reveals macro-level trends, while
the proportional map provides micro-level insights.
Input: Two selectInput() components allow users to:
(1) choose the year (2014–2024), and
(2) search and select a country for focused exploration.
Output: A coordinated view using two maps:
A choropleth map rendered via tmap, showing happiness levels by country with zoom-to-country feature.
A proportional symbol map rendered with leaflet, displaying happiness scores with interactive popups.
UI Components Used: selectInput() (for year and country), sidebarPanel() for inputs, and tmapOutput() + leafletOutput() within a fluidRow() layout for side-by-side map rendering.

# Load libraries
library(shiny); library(tidyverse); library(sf); library(tmap); library(leaflet); library(rnaturalearth)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year %in% 2014:2024)
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
left_join(happiness, by = c("name" = "country"))
ui <- fluidPage(
titlePanel("🌍 World Happiness Explorer"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:", choices = sort(unique(world_happy$year)), selected = 2024),
selectInput("selected_country", "Search Country:", choices = NULL)
),
mainPanel(
fluidRow(
column(6,
h4("Choropleth Map"),
tmapOutput("choropleth_map", height = "500px")
),
column(6,
h4("Proportional Symbol Map"),
leafletOutput("prop_map", height = "500px")
)
)
)
)
)
server <- function(input, output, session) {
# Reactive filtered data by year
filtered_data <- reactive({
world_happy %>%
filter(year == input$selected_year & !is.na(ladder_score))
})
# Update dropdown
observe({
updateSelectInput(session, "selected_country",
choices = sort(unique(filtered_data()$name)))
})
# Choropleth Map - Zoom to selected country
output$choropleth_map <- renderTmap({
tmap_mode("view")
# Zooming to selected country
selected_geom <- filtered_data() %>% filter(name == input$selected_country)
bbox_zoom <- if (nrow(selected_geom) > 0) st_bbox(selected_geom) else st_bbox(filtered_data())
tm_shape(filtered_data(), bbox = bbox_zoom) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name",
popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
title = paste("Happiness Score:", input$selected_year)
)
})
# Proportional Symbol Map
output$prop_map <- renderLeaflet({
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
leaflet(coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
})
# Auto-zoom on dropdown change
observeEvent(input$selected_country, {
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
selected_data <- coords %>% filter(name == input$selected_country)
leafletProxy("prop_map") %>%
setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
clearPopups() %>%
addPopups(
lng = selected_data$lon,
lat = selected_data$lat,
popup = paste0(
"<b>Country:</b> ", selected_data$name, "<br/>",
"<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
"<b>Region:</b> ", selected_data$region
)
)
})
# Auto-update dropdown when user clicks map
observeEvent(input$choropleth_map_shape_click, {
clicked_country <- input$choropleth_map_shape_click$id
updateSelectInput(session, "selected_country", selected = clicked_country)
})
}
shinyApp(ui, server)In this iteration, I introduced an additional selectInput() to filter the data by region. This allows users to explore happiness scores by broader geographical areas before narrowing down to a specific country. The country dropdown is dynamically updated based on the selected region, and selecting a country will also auto-update the region field—ensuring a smooth and linked filtering experience.
Input: Three selectInput() components allow users to:
(1) choose the year (2014–2024),
(2) filter by region,
(3) search and select a country within a selected region.
Output: A coordinated view using two maps:
A choropleth map rendered via tmap, showing happiness levels by country with zoom-to-country feature.
A proportional symbol map rendered with leaflet, displaying happiness scores with interactive popups.
UI Components Used: selectInput() (for year, region and country), sidebarPanel() for inputs, and tmapOutput() + leafletOutput() within a fluidRow() layout for side-by-side map rendering.

# Load libraries
library(shiny); library(tidyverse); library(sf); library(tmap); library(leaflet); library(rnaturalearth); library(readr)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year %in% 2014:2024)
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
left_join(happiness, by = c("name" = "country"))
# UI
ui <- fluidPage(
titlePanel("🌍 World Happiness Explorer"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:",
choices = sort(unique(world_happy$year)),
selected = 2024),
selectInput("selected_region", "Filter by Region:", choices = NULL),
selectInput("selected_country", "Search Country:", choices = NULL)
),
mainPanel(
fluidRow(
column(6,
h4("Choropleth Map"),
tmapOutput("choropleth_map", height = "500px")
),
column(6,
h4("Proportional Symbol Map"),
leafletOutput("prop_map", height = "500px")
)
)
)
)
)
# Server
server <- function(input, output, session) {
# Initial values for dropdowns
observe({
regions <- sort(unique(world_happy$region))
updateSelectInput(session, "selected_region",
choices = c("All", regions),
selected = "All")
})
# Reactive filtered data
filtered_data <- reactive({
data <- world_happy %>%
filter(year == input$selected_year & !is.na(ladder_score))
if (input$selected_region != "All") {
data <- data %>% filter(region == input$selected_region)
}
return(data)
})
# Update country dropdown when region changes
observeEvent(input$selected_region, {
countries <- filtered_data() %>% pull(name) %>% unique() %>% sort()
updateSelectInput(session, "selected_country",
choices = countries,
selected = countries[1])
})
# Update region when country changes
observeEvent(input$selected_country, {
selected_region <- world_happy %>%
filter(name == input$selected_country,
year == input$selected_year) %>%
pull(region) %>% unique()
if (!is.null(selected_region) && length(selected_region) == 1) {
updateSelectInput(session, "selected_region",
selected = selected_region)
}
})
# Choropleth Map
output$choropleth_map <- renderTmap({
tmap_mode("view")
selected_geom <- filtered_data() %>% filter(name == input$selected_country)
bbox_zoom <- if (nrow(selected_geom) > 0) st_bbox(selected_geom) else st_bbox(filtered_data())
tm_shape(filtered_data(), bbox = bbox_zoom) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name",
popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
title = paste("Happiness Score:", input$selected_year)
)
})
# Proportional Symbol Map
output$prop_map <- renderLeaflet({
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
leaflet(coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
})
# Auto-zoom and popup on country selection
observeEvent(input$selected_country, {
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
selected_data <- coords %>% filter(name == input$selected_country)
if (nrow(selected_data) > 0 &&
!is.na(selected_data$lon) &&
!is.na(selected_data$lat)) {
leafletProxy("prop_map") %>%
setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
clearPopups() %>%
addPopups(
lng = selected_data$lon,
lat = selected_data$lat,
popup = paste0(
"<b>Country:</b> ", selected_data$name, "<br/>",
"<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
"<b>Region:</b> ", selected_data$region
)
)
}
})
}
# Run the app
shinyApp(ui, server)In the next tab of our Cluster and Outlier Analysis module, I will apply spatial autocorrelation techniques to the World Happiness dataset by generating a Moran Scatter Plot and a Local Indicators of Spatial Association (LISA) Cluster Map.
The Moran Scatter Plot will quantify global spatial autocorrelation, indicating the overall degree of spatial clustering of happiness scores. The LISA Cluster Map will identify statistically significant local clusters—such as High-High (hot spots), Low-Low (cold spots), and spatial outliers (High-Low, Low-High)—allowing for granular analysis of localized spatial dependencies. T
The Moran scatter plot visualizes the spatial autocorrelation between countries’ Happiness Scores and the average scores of their neighboring countries. This helps us understand whether similar levels of happiness cluster geographically.
In the example below, we use moran.plot() from the spdep package to explore spatial relationships in the 2024 Happiness Score:
Interpretation:
The plot is split into four quadrants, each representing a different spatial relationship:
High-High (HH) - Top right quadrant: countries with high
Low-Low (LL) - Bottom-left quadrant: countries with low happiness surrounded by similarly unhappy neighbors.
High-Low (HL) - Bottom-right quadrant: — countries with high happiness surrounded by low-scoring neighbors (potential outliers).
Low-High (LH) - Top-left: countries with low happiness surrounded by happy neighbors (also outliers).
Note: Scatter plot provides an intuitive way to detect potential spatial clusters/ anomalies.
# Filter and prepare data
world_2024 <- world_happy %>% filter(year == 2024 & !is.na(ladder_score))
# Create spatial weights (Queen's case)
coords <- st_coordinates(st_centroid(world_2024))
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Moran Scatter Plot
moran.plot(world_2024$ladder_score, lw,
labels = as.character(world_2024$name),
xlab = "Happiness Score (2024)",
ylab = "Spatially Lagged Happiness Score",
zero.policy = TRUE)
Before generating the LISA cluster map, it is important to center both the spatially lagged variable and the local Moran’s I values. This step allows us to identify meaningful spatial patterns by determining whether a country’s happiness score — and the scores of its neighbors — are above or below the global mean.
By doing so, we can categorize each country into one of the four LISA cluster types (High-High, Low-Low, High-Low, Low-High). These classifications are only meaningful when assessed relative to the dataset’s overall distribution, making centering a crucial step for accurate spatial interpretation.
# Load required library
library(spdep)
# Step 1: Filter 2024 data
world_2024 <- world_happy %>%
filter(year == 2024 & !is.na(ladder_score))
# Step 2: Create spatial weights
coords <- st_coordinates(st_centroid(world_2024))
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Step 3: Calculate local Moran's I
localMI <- localmoran(world_2024$ladder_score, lw, zero.policy = TRUE)
# Step 4: Spatial lag of the variable
world_2024$lag_score <- lag.listw(lw, world_2024$ladder_score)
# Step 5: Center spatial lag & local Moran's I
DV <- world_2024$lag_score - mean(world_2024$lag_score, na.rm = TRUE)
LM_I <- localMI[, 1] - mean(localMI[, 1], na.rm = TRUE)
# Step 6: Set significance threshold
signif <- 0.05
# Step 7: Create quadrant classification
quadrant <- vector(mode = "numeric", length = nrow(localMI))
quadrant[DV < 0 & LM_I > 0] <- 1 # Low-Low
quadrant[DV > 0 & LM_I < 0] <- 2 # High-Low
quadrant[DV < 0 & LM_I < 0] <- 3 # Low-High
quadrant[DV > 0 & LM_I > 0] <- 4 # High-High
quadrant[localMI[, 5] > signif] <- 0 # Not significant
# Step 8: Assign quadrant to data
world_2024$quadrant <- quadrantThe LISA Cluster Map provides a powerful visualization of localized spatial autocorrelation in happiness scores for 2024. By computing the local Moran’s I for each country, we can identify areas where similar values cluster together (e.g., high-scoring countries surrounded by other high scorers), or where outliers emerge (e.g., a low-scoring country surrounded by high scorers).
Each country is classified into one of five categories based on the direction and significance of the spatial relationship:
| Category | Description | Statistically significant? | Interpretation |
|---|---|---|---|
| HH | High with high neighbours | Yes | Cluster (hot spot) |
| LL | Low with low neighbours | Yes | Cluster (cold spot) |
| HL | High with low neighbours | Yes | Spatial outlier |
| LH | Low with high neighbours | Yes | Spatial outlier |
| Insignificant | No strong autocorrelation | No | Insignifcant |
# Load libraries
library(tidyverse); library(sf); library(spdep); library(tmap); library(readr); library(rnaturalearth)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year == 2024 & !is.na(ladder_score))
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
# Join spatial and happiness data
world_2024 <- world %>%
left_join(happiness, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
# Construct spatial weights
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Compute local Moran's I
localMI <- localmoran(world_2024$ladder_score, lw, zero.policy = TRUE)
# Prepare LISA categories
quadrant <- vector(mode = "numeric", length = nrow(localMI))
lagged_score <- lag.listw(lw, world_2024$ladder_score)
centered_lag <- lagged_score - mean(lagged_score)
centered_localMI <- localMI[, 1] - mean(localMI[, 1])
significance_level <- 0.05
# Define cluster categories
quadrant[centered_lag < 0 & centered_localMI > 0] <- 1 # Low-Low
quadrant[centered_lag > 0 & centered_localMI < 0] <- 2 # Low-High
quadrant[centered_lag < 0 & centered_localMI < 0] <- 3 # High-Low
quadrant[centered_lag > 0 & centered_localMI > 0] <- 4 # High-High
quadrant[localMI[, 5] > significance_level] <- 0 # Not significant
# Add to spatial data
world_2024$quadrant <- quadrant
world_2024$cluster_label <- factor(
quadrant,
levels = 0:4,
labels = c("Insignificant", "Low-Low", "Low-High", "High-Low", "High-High")
)
# Plot the interactive LISA cluster map using tmap
tmap_mode("view")
tm_shape(world_2024) +
tm_fill(
col = "cluster_label",
palette = c(
"Insignificant" = "#ffffcc",
"Low-Low" = "blue",
"Low-High" = "#78c679",
"High-Low" = "#c2e699",
"High-High" = "red"
),
title = "LISA Cluster (2024)",
style = "cat",
id = "admin", # <<-- This makes popups/hover show country names
popup.vars = c(
"Country" = "admin",
"Cluster Type" = "cluster_label",
"Happiness Score" = "ladder_score"
)
) +
tm_borders(alpha = 0.4) +
tm_layout(
frame = FALSE,
legend.outside = TRUE
)# Load libraries
library(shiny); library(tidyverse); library(sf); library(spdep); library(tmap); library(readr); library(rnaturalearth)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(!is.na(ladder_score))
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
# UI
ui <- navbarPage("World Happiness Explorer",
tabPanel("Spatial Clustering",
fluidPage(
titlePanel("🌐 Spatial Clustering Analysis"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:",
choices = sort(unique(happiness$year)),
selected = 2024),
hr(),
h4("Chart Interpretation"),
helpText("The Moran scatterplot shows how each country's happiness score correlates with its neighbors'."),
helpText("The LISA Cluster map highlights statistically significant spatial clusters:",
"\n- High-High: Happy countries near other happy countries",
"\n- Low-Low: Unhappy countries near unhappy neighbors",
"\n- High-Low / Low-High: Potential outliers",
"\n- Insignificant: No strong spatial pattern")
),
mainPanel(
fluidRow(
column(6, plotOutput("moran_plot", height = "500px")),
column(6, tmapOutput("lisa_map", height = "500px"))
)
)
)
)
)
)
# Server
server <- function(input, output, session) {
world_data <- reactive({
data <- happiness %>%
filter(year == input$selected_year)
world %>%
left_join(data, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
})
# Spatial weights and local moran
local_moran <- reactive({
data <- world_data()
nb <- poly2nb(data, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
localMI <- localmoran(data$ladder_score, lw, zero.policy = TRUE)
lagged_score <- lag.listw(lw, data$ladder_score)
centered_lag <- lagged_score - mean(lagged_score)
centered_localMI <- localMI[, 1] - mean(localMI[, 1])
quadrant <- vector(mode = "numeric", length = nrow(data))
quadrant[centered_lag < 0 & centered_localMI > 0] <- 1 # Low-Low
quadrant[centered_lag > 0 & centered_localMI < 0] <- 2 # Low-High
quadrant[centered_lag < 0 & centered_localMI < 0] <- 3 # High-Low
quadrant[centered_lag > 0 & centered_localMI > 0] <- 4 # High-High
quadrant[localMI[, 5] > 0.05] <- 0 # Not significant
data$quadrant <- quadrant
data$cluster_label <- factor(
quadrant,
levels = 0:4,
labels = c("Insignificant", "Low-Low", "Low-High", "High-Low", "High-High")
)
list(data = data, lw = lw)
})
# Moran plot
output$moran_plot <- renderPlot({
dat <- local_moran()
moran.plot(dat$data$ladder_score, dat$lw,
labels = dat$data$name,
xlab = paste("Happiness Score (", input$selected_year, ")", sep = ""),
ylab = "Spatially Lagged Happiness Score",
zero.policy = TRUE)
})
# LISA map
output$lisa_map <- renderTmap({
tmap_mode("view")
tm_shape(local_moran()$data) +
tm_fill(
col = "cluster_label",
palette = c(
"Insignificant" = "#ffffcc",
"Low-Low" = "blue",
"Low-High" = "#78c679",
"High-Low" = "#c2e699",
"High-High" = "red"
),
title = paste("LISA Cluster (", input$selected_year, ")", sep = ""),
style = "cat",
id = "admin",
popup.vars = c(
"Country" = "admin",
"Cluster Type" = "cluster_label",
"Happiness Score" = "ladder_score"
)
) +
tm_borders(alpha = 0.4) +
tm_layout(frame = FALSE, legend.outside = TRUE)
})
}
shinyApp(ui = ui, server = server)---
title: "Take-home Exercise 03"
author: "Andrea Yeo"
date-modified: "last-modified"
date: "March 22, 2025"
execute:
echo: true
eval: true
warning: false
freeze: true
---
# Prototyping Modules for Visual Analytics Shiny Application
## 1. The task
For this take-home exercise, each team member have to choose one module from our proposed Shiny application and complete the following tasks:
1. **Verify package support:** Identify and ensure all required R packages for the module are available on CRAN.
2. **Code testing**: Prepare and test the relevant R code to confirm it runs correctly and produces the expected output.
3. **Define inputs and outputs**: Specify the parameters (inputs) and outputs that will be featured in the Shiny application.
4. **Select UI components**: Choose suitable Shiny UI elements to present and control the identified parameters within the app interface.
## 2. Getting started
Our project utilizes open-source data from the [World Happiness Report – Data Sharing](https://worldhappiness.report/data-sharing/) covering the years **2015 to 2025**. Detailed information on the data cleaning and preparation process is available on our group's Netlify site, under the **"Data Preparation"** tab: [Group Visual Analytic Application](https://worldhappiness-vaa.netlify.app/).
Specifically, our project focuses on the visualisation of global happiness trends. I will be exploring the modules related to **geospatial** and **aspatial analysis** to explore both spatial patterns and non-spatial factors influencing happiness scores across countries. For the prototyping of these modules, I will be using the **`cleaned_data`** dataset prepared after the data processing phase.
### 2.1 Loading R packages
The below R packages will be used in this exercise and for the Shiny application
```{r}
pacman::p_load(sf, tidyverse, tmap, dplyr,
spatstat, spdep,
lubridate, leaflet,
plotly, DT, viridis,
ggplot2, sfdep, janitor, shiny, rnaturalearth)
```
### 2.2 Reading the dataset
```{r}
happiness <- read.csv("data/world_happiness.csv")
```
### 2.3 Understanding the data structure
- `glimpse()`: provides a transposed overview of a dataset, showing variables and their types in a concise format.
- `head()`: displays the first few rows of a dataset (default is 6 rows) to give a quick preview of the data.
- `summary()`: generates a statistical summary of each variable, including measures like mean, median, and range for numeric data.
- `duplicated()`:returns a logical vector indicating which elements or rows in a vector or data frame are duplicates.
- `colSums(is.na())`: counts the number of missing values (NA) in each column of the data frame.
- `str()`: use `str()` to display the column names, data types, and a preview of the data.
:::::: panel-tabset
## glimpse()
```{r}
glimpse(happiness)
```
## head()
```{r}
head(happiness)
```
## summary()
```{r}
summary(happiness)
```
## duplicated()
```{r}
happiness[duplicated(happiness),]
```
::: callout-note
- Ensure that there are no duplicated columns, if not will have to investigate further.
:::
## colSum(is.na())
```{r}
colSums(is.na(happiness))
```
::: callout-note
- Ensure that there are no NA values, if not will have to investigate further.
- Possibility to use `drop_na()` function to drop rows where any specified column contains a missing value.
:::
## str())
```{r}
str(happiness)
```
::: callout-note
- Ensure that all variables are correctly classified by data type; recast variable types if needed.
- Variables are correctly classified - where categorical variables are classified as **character**, while continuous variables are classified as **double**.
:::
::::::
The `happiness` tibble contains 12 attributes, as shown above.
The following preprocessing checks were conducted as part of data preparation:
::: callout-tip
## Preprocessing Checks
- Verified that the correct data types were loaded in the `happiness` dataset using `glimpse()` and `str()`
- Ensured there were no duplicate variable names using `duplicated()` in the dataset
- Checked for missing values using `colSums(is.na())`
:::
::: callout-tip
## Data scale review
- The column values are already **normalized**, ranging consistently between **1 and 10** across countries.
- No further standardization/ transformation is required, as the scale is uniform and interpretable for spatial analysis.
:::
## 3. Geospatial analysis (Choropleth Map)
### 3.1 To load world country boundaries
To support the geospatial visualisation of global happiness trends, the [`rnaturalearth`](https://cran.r-project.org/web/packages/rnaturalearth/vignettes/rnaturalearth.html) package was used to import spatial boundary data of countries. Specifically, the `ne_countries()` function was called with a medium scale and returned as a simple features (`sf`) object, stored in the `world` variable. This spatial dataset provides the geographic outlines of all countries and is essential for merging with the **World Happiness Report data** to enable choropleth mapping and spatial analysis at the country level.
```{r}
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
```
### 3.2 Join happiness data with geospatial features
To prepare the dataset for mapping, the World Happiness data was filtered to include only records from the year 2024. This filtered data (`happiness_latest`) was then joined with the spatial dataset (`world`) using the country name as the common key. The resulting dataset, `world_happy`, combines both geographic boundaries and happiness scores, enabling geospatial visualisation of happiness levels across countries.
```{r}
# Join by country name
happiness_latest <- happiness %>%
filter(year == 2024)
world_happy <- world %>%
left_join(happiness_latest, by = c("name" = "country"))
```
### 3.3 Choropleth map of global happiness (basic version)
The code below plots a choropleth map to visualise the distribution of happiness scores across countries for the year 2024. Using the [`tmap`](https://cran.r-project.org/web/packages/tmap/index.html) package in plotting mode, the `tm_shape()` function defines the spatial data object (`world_happy`), which contains both geographic boundaries and the associated happiness data. The `tm_polygons()` function maps the `ladder_score` variable to a blue color gradient, with darker shades indicating higher happiness scores.
Finally, `tm_layout()` adds a descriptive title and positions the legend outside the map area to improve clarity and presentation.
```{r}
tmap_mode("plot")
tmap_style("white")
tm_shape(world_happy) +
tm_polygons("ladder_score", palette = "Blues", title = "Happiness Score (2024)") +
tm_layout(title = "World Happiness Map", legend.outside = TRUE)
```
### 3.4 Choropleth map with styling and classification (enhanced version)
The enhanced choropleth map uses the [`tmap`](https://cran.r-project.org/web/packages/tmap/index.html) package with improved styling and classification. It applies a quantile-based scheme with five classes and the `"YlOrRd"` palette to represent `ladder_score`, where darker shades indicate higher happiness. Additional cartographic elements such as a compass, scale bar, borders, and descriptive credits are included to improve readability and interpretability. The legend is placed on the right, and a classic map style is applied for a cleaner layout.
```{r}
tmap_style("classic")
tmap_mode("plot")
tm_shape(world_happy) +
tm_polygons(
col = "ladder_score",
palette = "YlOrRd",
style = "quantile",
n = 5,
title = "Happiness Score (2024)"
) +
tm_layout(
main.title = "World Happiness Distribution (2024)",
main.title.position = "center",
main.title.size = 1.2,
legend.outside = TRUE,
legend.position = c("right", "center"),
frame = FALSE
) +
tm_borders(alpha = 0.4, lwd = 0.3) +
tm_compass(type = "8star", size = 2, position = c("left", "top")) +
tm_scale_bar(position = c("left", "bottom")) +
tm_credits("Source: World Happiness Report 2024", position = c("left", "bottom")) +
tm_credits("Note: Darker shades indicate higher happiness; lighter shades indicate lower happiness.",
position = c("left", "bottom"), just = "right", size = 0.6)
```
### 3.5 Summary of differences between Basic version and Enhanced Version
| Feature | Basic Version | Enhanced Version |
|------------------------|------------------------|------------------------|
| Classification method | Default("pretty") | Custom("quantile", 5 bins) |
| Layout | Basic title, default layout | Custom title, centered, larger size |
| Map Borders | None | Added with transparency and thin width |
| Compass/ Scale Bar | Not included | Included |
| Credits/ Source | Not included | Included |
| Frame | Default | Removed |
## 4. Drawing small multiple (faceted) map
The code below creates a faceted choropleth map to visualise global happiness scores from 2014 to 2024. Using [tmap](https://cran.r-project.org/web/packages/tmap/index.html), each facet represents a different year, allowing for easy comparison of happiness trends across countries. Geospatial data is joined using [rnaturalearth](https://cran.r-project.org/web/packages/rnaturalearth/vignettes/rnaturalearth.html).
```{r fig.width=16, fig.height=20}
# Load required libraries
library(tidyverse); library(tmap); library(sf); library(rnaturalearth)
# Set style and plotting mode
tmap_style("white")
tmap_mode("plot")
# Step 1: Filter for all years from 2014 to 2024
years_to_plot <- 2014:2024
happiness_subset <- happiness %>%
filter(year %in% years_to_plot)
# Step 2: Load and join with world spatial data
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy_facet <- world %>%
left_join(happiness_subset, by = c("name" = "country"))
# Step 3: Filter out rows with missing ladder_score
world_happy_facet_clean <- world_happy_facet %>%
filter(!is.na(ladder_score))
# Step 4: Plot small multiple choropleth maps
tm_shape(world_happy_facet_clean) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
style = "quantile",
title = "Happiness Score"
) +
tm_facets(by = "year", ncol = 2) +
tm_layout(
main.title = "Global Happiness Scores (2014–2024)",
main.title.position = "center",
panel.label.size = 2.0,
legend.outside = TRUE,
legend.position = c("center", "top"),
legend.height = 1.5,
legend.width = 1.0,
legend.text.size = 1.5,
legend.title.size = 1.5,
asp = 0,
between.margin = 0,
outer.margins = c(0, 0, 0, 0)
)
```
## 5. Animating global happiness over time
The code below produces an animated choropleth map using `ggplot2`, `gganimate`, and `sf` to illustrate the temporal evolution of global happiness scores from 2014 to 2024. Country-level geospatial data from `rnaturalearth` is merged with World Happiness Report data via a spatial join. `transition_time()` is used to animate by year, with countries shaded according to their ladder score.
::: callout-note
**Note:** Although the underlying world geometries are static, the animated map may appear to shift or jitter between frames. This occurs because some countries are missing data in certain years, leading to variations in the number of plotted geometries. As a result, `gganimate` recalculates the plot layout for each frame, even when `coord_sf(xlim, ylim, expand = FALSE)` is specified. The visual movement is not due to changes in the map itself, but rather to inconsistencies in data availability across years.
:::
```{r}
library(tidyverse); library(sf); library(gganimate); library(rnaturalearth)
# Load base world map
world <- ne_countries(scale = "medium", returnclass = "sf") %>%
select(name, geometry)
# Load and clean happiness data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year %in% 2014:2024)
# Merge first — keep all geometries
world_happy <- world %>%
left_join(happiness, by = c("name" = "country")) %>%
mutate(year = as.numeric(year))
# Plot
ggplot(world_happy) +
geom_sf(aes(fill = ladder_score)) +
scale_fill_distiller(palette = "YlGnBu", direction = 1, na.value = "lightgray") +
labs(title = "Happiness Score: {frame_time}", fill = "Score") +
transition_time(year) +
ease_aes('linear') +
coord_sf(
xlim = c(-180, 180),
ylim = c(-60, 90),
expand = FALSE
)
```
## 6. Exploring global happiness with an interactive map - Choropleth Map
This code creates an interactive world map using `tmap`, showing happiness scores by country for 2024. Each country is colored based on its ladder score, and tooltips display the country name when hovered. The map allows users to zoom and explore happiness data visually.
```{r}
happiness_latest <- happiness_subset %>%
filter(year == 2024)
world_happy_latest <- world %>%
left_join(happiness_latest, by = c("name" = "country"))
world_happy_latest_clean <- world_happy_latest %>%
filter(!is.na(ladder_score))
tmap_mode("view")
tm_shape(world_happy_latest_clean) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name", # this is what shows in the tooltip
title = "Happiness Score (2024)"
)
```
## 7. User interface design - Part 1 - Choropleth Map
In this module, I developed an interactive Shiny application to visualize the World Happiness Score by country and year. The key **input** is a dropdown menu `(selectInput)` that allows users to choose a specific year from 2014 to 2024. Based on the selected year, the **output** is a dynamic choropleth map `(tmapOutput)` showing each country's happiness score using a color gradient.
::: callout-note
- **Input:** A single `selectInput()` allows the user to choose the **Year** of interest from 2014 to 2024.
- **Output:** A dynamic choropleth map is rendered using `tmap`, updating automatically based on the selected year.
- **UI Components Used:** `selectInput()` for the dropdown menu, `sidebarPanel()` for input placement, and `tmapOutput()` in `mainPanel()` for rendering the interactive map.
:::
::: panel-tabset
## UI()

## Code()
```{r}
# Load libraries
library(shiny); library(tmap); library(tidyverse); library(sf); library(rnaturalearth)
# Prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year))
years_to_plot <- 2014:2024
happiness_subset <- happiness %>%
filter(year %in% years_to_plot)
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
left_join(happiness_subset, by = c("name" = "country"))
# UI
ui <- fluidPage(
titlePanel("🌍 World Happiness Score by Year"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:",
choices = sort(unique(world_happy$year)),
selected = 2024)
),
mainPanel(
tmapOutput("happiness_map")
)
)
)
# Server
server <- function(input, output, session) {
filtered_data <- reactive({
world_happy %>%
filter(year == input$selected_year & !is.na(ladder_score))
})
output$happiness_map <- renderTmap({
tmap_mode("view")
tm_shape(filtered_data()) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name",
popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
title = paste("Happiness Score:", input$selected_year)
)
})
}
# Run the app
shinyApp(ui, server)
```
:::
## 8. Mapping global happiness with proportional circles - Proportional Symbol Map
This code below generates an interactive map using `leaflet` to visualize 2024 World Happiness Scores. Each country is displayed as a circle located at its geographic center, with the size representing its happiness score. I added tooltips which provide additional insights like economy, life expectancy, freedom, and region, and a legend which helps interpret the happiness scores visually.
```{r}
library(tidyverse); library(readr); library(leaflet)
library(sf); library(rnaturalearth); library(rnaturalearthdata)
# Load the world happiness data
happiness <- read_csv("data/world_happiness.csv")
# Filter for a specific year (e.g., 2024)
happiness_latest <- happiness %>%
filter(year == 2024)
# Load country polygons and calculate centroids
world_sf <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
select(name, geometry)
# Join happiness data to spatial data
world_happy_sf <- world_sf %>%
left_join(happiness_latest, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
# Get centroids for each country
world_happy_points <- st_centroid(world_happy_sf)
# Extract coordinates for leaflet
world_happy_points_coords <- world_happy_points %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
leaflet(world_happy_points_coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
```
## 9. User interface design - Part 2 - Proportional Symbol Map
In this module, I developed an interactive Shiny application that allows users to explore the World Happiness Score by country for the year 2024. The main **input** is a country selector using `selectInput`, which enables users to search and select a country from the dropdown list. Once selected, the **output** is an interactive `leaflet` map that **automatically zooms** to the selected country and displays a popup with additional details such as economy, life expectancy, freedom, and region.
::: callout-note
- **Input:** A `selectInput()` enables users to search and select a country from the list of available options for 2024.
- **Output:** A proportional symbol map rendered via `leaflet`, displaying country-level happiness scores using interactive circle markers with popups.
- **UI Components Used:** `selectInput()` in the `sidebarPanel()` for country selection, and `leafletOutput()` in the `mainPanel()` for displaying the dynamic map.
:::
::: panel-tabset
## UI()

## Code()
```{r}
library(shiny); library(tidyverse); library(readr); library(leaflet)
library(sf); library(rnaturalearth); library(rnaturalearthdata)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
filter(year == 2024)
world_sf <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
select(name, geometry)
world_happy_sf <- world_sf %>%
left_join(happiness, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
world_happy_points <- st_centroid(world_happy_sf)
world_happy_points_coords <- world_happy_points %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
# UI
ui <- fluidPage(
titlePanel("World Happiness (2024)"),
sidebarLayout(
sidebarPanel(
selectInput("selected_country", "Search Country:",
choices = sort(world_happy_points_coords$name),
selected = NULL)
),
mainPanel(
leafletOutput("map", height = "400px")
)
)
)
# Server
server <- function(input, output, session) {
output$map <- renderLeaflet({
leaflet(world_happy_points_coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
})
observeEvent(input$selected_country, {
selected_data <- world_happy_points_coords %>%
filter(name == input$selected_country)
leafletProxy("map") %>%
setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
clearPopups() %>%
addPopups(
lng = selected_data$lon,
lat = selected_data$lat,
popup = paste0(
"<b>Country:</b> ", selected_data$name, "<br/>",
"<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
"<b>Region:</b> ", selected_data$region
)
)
})
}
# Run the app
shinyApp(ui, server)
```
:::
## 10. User interface design - Part 3 - Combining both Choropleth Map & Proportional Symbol Map
This module presents an interactive Shiny application with two coordinated maps to explore World Happiness data. The **Choropleth Map** allows users to select a year and view global happiness patterns using a color gradient. In contrast, the **Proportional Symbol Map** lets users search for a country, zoom in automatically, and view detailed indicators such as economy, life expectancy, and freedom. Circle size and color reflect the happiness score.
After considering various approaches, I chose to implement both map types as they serve different purposes:
- the choropleth reveals macro-level trends, while
- the proportional map provides micro-level insights.
::: callout-tip
- Note: This application includes smart interactivity features. Selecting a country from the dropdown automatically **zooms both maps** to the selected country for focused viewing. Conversely, **clicking a country** on the choropleth map will **update the dropdown menu** —ensuring two-way interaction. These features improve usability by making geographic exploration more intuitive and dynamic.
:::
::: callout-note
- **Input:** Two `selectInput()` components allow users to:
- \(1\) choose the year (2014–2024), and
- \(2\) search and select a country for focused exploration.
- **Output:** A *coordinated view* using two maps:
- A choropleth map rendered via `tmap,` showing happiness levels by country with zoom-to-country feature.
- A proportional symbol map rendered with `leaflet`, displaying happiness scores with interactive popups.
- **UI Components Used:** `selectInput()` (for year and country), `sidebarPanel()` for inputs, and `tmapOutput()` + `leafletOutput()` within a `fluidRow()` layout for side-by-side map rendering.
:::
::: panel-tabset
## UI()

## Code()
```{r}
# Load libraries
library(shiny); library(tidyverse); library(sf); library(tmap); library(leaflet); library(rnaturalearth)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year %in% 2014:2024)
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
left_join(happiness, by = c("name" = "country"))
ui <- fluidPage(
titlePanel("🌍 World Happiness Explorer"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:", choices = sort(unique(world_happy$year)), selected = 2024),
selectInput("selected_country", "Search Country:", choices = NULL)
),
mainPanel(
fluidRow(
column(6,
h4("Choropleth Map"),
tmapOutput("choropleth_map", height = "500px")
),
column(6,
h4("Proportional Symbol Map"),
leafletOutput("prop_map", height = "500px")
)
)
)
)
)
server <- function(input, output, session) {
# Reactive filtered data by year
filtered_data <- reactive({
world_happy %>%
filter(year == input$selected_year & !is.na(ladder_score))
})
# Update dropdown
observe({
updateSelectInput(session, "selected_country",
choices = sort(unique(filtered_data()$name)))
})
# Choropleth Map - Zoom to selected country
output$choropleth_map <- renderTmap({
tmap_mode("view")
# Zooming to selected country
selected_geom <- filtered_data() %>% filter(name == input$selected_country)
bbox_zoom <- if (nrow(selected_geom) > 0) st_bbox(selected_geom) else st_bbox(filtered_data())
tm_shape(filtered_data(), bbox = bbox_zoom) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name",
popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
title = paste("Happiness Score:", input$selected_year)
)
})
# Proportional Symbol Map
output$prop_map <- renderLeaflet({
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
leaflet(coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
})
# Auto-zoom on dropdown change
observeEvent(input$selected_country, {
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
selected_data <- coords %>% filter(name == input$selected_country)
leafletProxy("prop_map") %>%
setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
clearPopups() %>%
addPopups(
lng = selected_data$lon,
lat = selected_data$lat,
popup = paste0(
"<b>Country:</b> ", selected_data$name, "<br/>",
"<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
"<b>Region:</b> ", selected_data$region
)
)
})
# Auto-update dropdown when user clicks map
observeEvent(input$choropleth_map_shape_click, {
clicked_country <- input$choropleth_map_shape_click$id
updateSelectInput(session, "selected_country", selected = clicked_country)
})
}
shinyApp(ui, server)
```
:::
## 11. User interface design - Part 4 - Adding regions
In this iteration, I introduced an additional `selectInput()` to filter the data by **region**. This allows users to explore happiness scores by broader geographical areas before narrowing down to a specific country. The country dropdown is dynamically updated based on the selected region, and selecting a country will also auto-update the region field—ensuring a smooth and linked filtering experience.
::: callout-tip
- Note: To explore and ensure that the *region* and *country* drop downs are dynamically linked, where
- selecting a *region* will show only countries in that region
- selecting a *country* will auto-update the region selection to match
:::
::: callout-note
- **Input:** Three `selectInput()` components allow users to:
- \(1\) choose the year (2014–2024),
- \(2\) filter by region,
- \(3\) search and select a country within a selected region.
- **Output:** A *coordinated view* using two maps:
- A choropleth map rendered via `tmap,` showing happiness levels by country with zoom-to-country feature.
- A proportional symbol map rendered with `leaflet`, displaying happiness scores with interactive popups.
- **UI Components Used:** `selectInput()` (for year, region and country), `sidebarPanel()` for inputs, and `tmapOutput()` + `leafletOutput()` within a `fluidRow()` layout for side-by-side map rendering.
:::
::: panel-tabset
## UI()

## Code()
```{r}
# Load libraries
library(shiny); library(tidyverse); library(sf); library(tmap); library(leaflet); library(rnaturalearth); library(readr)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year %in% 2014:2024)
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
left_join(happiness, by = c("name" = "country"))
# UI
ui <- fluidPage(
titlePanel("🌍 World Happiness Explorer"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:",
choices = sort(unique(world_happy$year)),
selected = 2024),
selectInput("selected_region", "Filter by Region:", choices = NULL),
selectInput("selected_country", "Search Country:", choices = NULL)
),
mainPanel(
fluidRow(
column(6,
h4("Choropleth Map"),
tmapOutput("choropleth_map", height = "500px")
),
column(6,
h4("Proportional Symbol Map"),
leafletOutput("prop_map", height = "500px")
)
)
)
)
)
# Server
server <- function(input, output, session) {
# Initial values for dropdowns
observe({
regions <- sort(unique(world_happy$region))
updateSelectInput(session, "selected_region",
choices = c("All", regions),
selected = "All")
})
# Reactive filtered data
filtered_data <- reactive({
data <- world_happy %>%
filter(year == input$selected_year & !is.na(ladder_score))
if (input$selected_region != "All") {
data <- data %>% filter(region == input$selected_region)
}
return(data)
})
# Update country dropdown when region changes
observeEvent(input$selected_region, {
countries <- filtered_data() %>% pull(name) %>% unique() %>% sort()
updateSelectInput(session, "selected_country",
choices = countries,
selected = countries[1])
})
# Update region when country changes
observeEvent(input$selected_country, {
selected_region <- world_happy %>%
filter(name == input$selected_country,
year == input$selected_year) %>%
pull(region) %>% unique()
if (!is.null(selected_region) && length(selected_region) == 1) {
updateSelectInput(session, "selected_region",
selected = selected_region)
}
})
# Choropleth Map
output$choropleth_map <- renderTmap({
tmap_mode("view")
selected_geom <- filtered_data() %>% filter(name == input$selected_country)
bbox_zoom <- if (nrow(selected_geom) > 0) st_bbox(selected_geom) else st_bbox(filtered_data())
tm_shape(filtered_data(), bbox = bbox_zoom) +
tm_polygons(
col = "ladder_score",
palette = "YlGnBu",
id = "name",
popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
title = paste("Happiness Score:", input$selected_year)
)
})
# Proportional Symbol Map
output$prop_map <- renderLeaflet({
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
leaflet(coords) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addCircleMarkers(
lng = ~lon,
lat = ~lat,
radius = ~ladder_score * 3,
color = "black",
fillColor = ~colorNumeric("YlGnBu", domain = coords$ladder_score)(ladder_score),
fillOpacity = 0.6,
stroke = TRUE,
weight = 0.5,
popup = ~paste0(
"<b>Country:</b> ", name, "<br/>",
"<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
"<b>Region:</b> ", region
)
) %>%
addLegend(
"bottomright",
pal = colorNumeric("YlGnBu", domain = coords$ladder_score),
values = ~ladder_score,
title = "Happiness Score",
opacity = 1
)
})
# Auto-zoom and popup on country selection
observeEvent(input$selected_country, {
centroids <- st_centroid(filtered_data())
coords <- centroids %>%
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
)
selected_data <- coords %>% filter(name == input$selected_country)
if (nrow(selected_data) > 0 &&
!is.na(selected_data$lon) &&
!is.na(selected_data$lat)) {
leafletProxy("prop_map") %>%
setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
clearPopups() %>%
addPopups(
lng = selected_data$lon,
lat = selected_data$lat,
popup = paste0(
"<b>Country:</b> ", selected_data$name, "<br/>",
"<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
"<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
"<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
"<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
"<b>Region:</b> ", selected_data$region
)
)
}
})
}
# Run the app
shinyApp(ui, server)
```
:::
## 12. Exploring the LISA cluster map
In the next tab of our *Cluster and Outlier Analysis* module, I will apply spatial autocorrelation techniques to the World Happiness dataset by generating a Moran Scatter Plot and a Local Indicators of Spatial Association (LISA) Cluster Map.
The Moran Scatter Plot will quantify global spatial autocorrelation, indicating the overall degree of spatial clustering of happiness scores. The LISA Cluster Map will identify statistically significant local clusters—such as High-High (hot spots), Low-Low (cold spots), and spatial outliers (High-Low, Low-High)—allowing for granular analysis of localized spatial dependencies. T
### 12.1 Plotting the Moran Scatter Plot
The **Moran scatter plot** visualizes the spatial autocorrelation between countries' *Happiness Scores* and the average scores of their neighboring countries. This helps us understand whether similar levels of happiness cluster geographically.
In the example below, we use `moran.plot()` from the [**spdep**](https://cran.r-project.org/web/packages/spdep/spdep.pdf) package to explore spatial relationships in the 2024 Happiness Score:
**Interpretation:**
The plot is split into **four quadrants**, each representing a different spatial relationship:
- **High-High (HH) - Top right quadrant:** countries with high
- **Low-Low (LL) -** **Bottom-left quadrant:** countries with low happiness surrounded by similarly unhappy neighbors.
- **High-Low (HL) - Bottom-right quadrant:** — countries with high happiness surrounded by low-scoring neighbors (potential outliers).
- **Low-High (LH) - Top-left:** countries with low happiness surrounded by happy neighbors (also outliers).
::: callout-tip
Note: Scatter plot provides an intuitive way to detect potential spatial clusters/ anomalies.
:::
```{r}
# Filter and prepare data
world_2024 <- world_happy %>% filter(year == 2024 & !is.na(ladder_score))
# Create spatial weights (Queen's case)
coords <- st_coordinates(st_centroid(world_2024))
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Moran Scatter Plot
moran.plot(world_2024$ladder_score, lw,
labels = as.character(world_2024$name),
xlab = "Happiness Score (2024)",
ylab = "Spatially Lagged Happiness Score",
zero.policy = TRUE)
```
### 12.2 Preparing LISA cluster class for Word Happiness (2024)
Before generating the LISA cluster map, it is important to center both the spatially lagged variable and the local Moran’s I values. This step allows us to identify meaningful spatial patterns by determining whether a country’s happiness score — and the scores of its neighbors — are above or below the global mean.
By doing so, we can categorize each country into one of the four LISA cluster types **(High-High, Low-Low, High-Low, Low-High)**. These classifications are only meaningful when assessed relative to the dataset’s overall distribution, making centering a crucial step for accurate spatial interpretation.
```{r}
# Load required library
library(spdep)
# Step 1: Filter 2024 data
world_2024 <- world_happy %>%
filter(year == 2024 & !is.na(ladder_score))
# Step 2: Create spatial weights
coords <- st_coordinates(st_centroid(world_2024))
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Step 3: Calculate local Moran's I
localMI <- localmoran(world_2024$ladder_score, lw, zero.policy = TRUE)
# Step 4: Spatial lag of the variable
world_2024$lag_score <- lag.listw(lw, world_2024$ladder_score)
# Step 5: Center spatial lag & local Moran's I
DV <- world_2024$lag_score - mean(world_2024$lag_score, na.rm = TRUE)
LM_I <- localMI[, 1] - mean(localMI[, 1], na.rm = TRUE)
# Step 6: Set significance threshold
signif <- 0.05
# Step 7: Create quadrant classification
quadrant <- vector(mode = "numeric", length = nrow(localMI))
quadrant[DV < 0 & LM_I > 0] <- 1 # Low-Low
quadrant[DV > 0 & LM_I < 0] <- 2 # High-Low
quadrant[DV < 0 & LM_I < 0] <- 3 # Low-High
quadrant[DV > 0 & LM_I > 0] <- 4 # High-High
quadrant[localMI[, 5] > signif] <- 0 # Not significant
# Step 8: Assign quadrant to data
world_2024$quadrant <- quadrant
```
### 12.3 Plotting the LISA Cluster Map
The LISA Cluster Map provides a powerful visualization of localized spatial autocorrelation in happiness scores for 2024. By computing the *local Moran’s I* for each country, we can identify areas where similar values cluster together (e.g., high-scoring countries surrounded by other high scorers), or where outliers emerge (e.g., a low-scoring country surrounded by high scorers).
Each country is classified into one of five categories based on the direction and significance of the spatial relationship:
| Category | Description | Statistically significant? | Interpretation |
|-------------|-----------------------|--------------|-----------------------|
| HH | High with high neighbours | Yes | Cluster (hot spot) |
| LL | Low with low neighbours | Yes | Cluster (cold spot) |
| HL | High with low neighbours | Yes | Spatial outlier |
| LH | Low with high neighbours | Yes | Spatial outlier |
| Insignificant | No strong autocorrelation | No | Insignifcant |
```{r}
# Load libraries
library(tidyverse); library(sf); library(spdep); library(tmap); library(readr); library(rnaturalearth)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(year == 2024 & !is.na(ladder_score))
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
# Join spatial and happiness data
world_2024 <- world %>%
left_join(happiness, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
# Construct spatial weights
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Compute local Moran's I
localMI <- localmoran(world_2024$ladder_score, lw, zero.policy = TRUE)
# Prepare LISA categories
quadrant <- vector(mode = "numeric", length = nrow(localMI))
lagged_score <- lag.listw(lw, world_2024$ladder_score)
centered_lag <- lagged_score - mean(lagged_score)
centered_localMI <- localMI[, 1] - mean(localMI[, 1])
significance_level <- 0.05
# Define cluster categories
quadrant[centered_lag < 0 & centered_localMI > 0] <- 1 # Low-Low
quadrant[centered_lag > 0 & centered_localMI < 0] <- 2 # Low-High
quadrant[centered_lag < 0 & centered_localMI < 0] <- 3 # High-Low
quadrant[centered_lag > 0 & centered_localMI > 0] <- 4 # High-High
quadrant[localMI[, 5] > significance_level] <- 0 # Not significant
# Add to spatial data
world_2024$quadrant <- quadrant
world_2024$cluster_label <- factor(
quadrant,
levels = 0:4,
labels = c("Insignificant", "Low-Low", "Low-High", "High-Low", "High-High")
)
# Plot the interactive LISA cluster map using tmap
tmap_mode("view")
tm_shape(world_2024) +
tm_fill(
col = "cluster_label",
palette = c(
"Insignificant" = "#ffffcc",
"Low-Low" = "blue",
"Low-High" = "#78c679",
"High-Low" = "#c2e699",
"High-High" = "red"
),
title = "LISA Cluster (2024)",
style = "cat",
id = "admin", # <<-- This makes popups/hover show country names
popup.vars = c(
"Country" = "admin",
"Cluster Type" = "cluster_label",
"Happiness Score" = "ladder_score"
)
) +
tm_borders(alpha = 0.4) +
tm_layout(
frame = FALSE,
legend.outside = TRUE
)
```
### 12.4 User interface design - Combining Moran Scattterplot with LISA map
```{r}
# Load libraries
library(shiny); library(tidyverse); library(sf); library(spdep); library(tmap); library(readr); library(rnaturalearth)
# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
mutate(year = as.numeric(year)) %>%
filter(!is.na(ladder_score))
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
# UI
ui <- navbarPage("World Happiness Explorer",
tabPanel("Spatial Clustering",
fluidPage(
titlePanel("🌐 Spatial Clustering Analysis"),
sidebarLayout(
sidebarPanel(
selectInput("selected_year", "Select Year:",
choices = sort(unique(happiness$year)),
selected = 2024),
hr(),
h4("Chart Interpretation"),
helpText("The Moran scatterplot shows how each country's happiness score correlates with its neighbors'."),
helpText("The LISA Cluster map highlights statistically significant spatial clusters:",
"\n- High-High: Happy countries near other happy countries",
"\n- Low-Low: Unhappy countries near unhappy neighbors",
"\n- High-Low / Low-High: Potential outliers",
"\n- Insignificant: No strong spatial pattern")
),
mainPanel(
fluidRow(
column(6, plotOutput("moran_plot", height = "500px")),
column(6, tmapOutput("lisa_map", height = "500px"))
)
)
)
)
)
)
# Server
server <- function(input, output, session) {
world_data <- reactive({
data <- happiness %>%
filter(year == input$selected_year)
world %>%
left_join(data, by = c("name" = "country")) %>%
filter(!is.na(ladder_score))
})
# Spatial weights and local moran
local_moran <- reactive({
data <- world_data()
nb <- poly2nb(data, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
localMI <- localmoran(data$ladder_score, lw, zero.policy = TRUE)
lagged_score <- lag.listw(lw, data$ladder_score)
centered_lag <- lagged_score - mean(lagged_score)
centered_localMI <- localMI[, 1] - mean(localMI[, 1])
quadrant <- vector(mode = "numeric", length = nrow(data))
quadrant[centered_lag < 0 & centered_localMI > 0] <- 1 # Low-Low
quadrant[centered_lag > 0 & centered_localMI < 0] <- 2 # Low-High
quadrant[centered_lag < 0 & centered_localMI < 0] <- 3 # High-Low
quadrant[centered_lag > 0 & centered_localMI > 0] <- 4 # High-High
quadrant[localMI[, 5] > 0.05] <- 0 # Not significant
data$quadrant <- quadrant
data$cluster_label <- factor(
quadrant,
levels = 0:4,
labels = c("Insignificant", "Low-Low", "Low-High", "High-Low", "High-High")
)
list(data = data, lw = lw)
})
# Moran plot
output$moran_plot <- renderPlot({
dat <- local_moran()
moran.plot(dat$data$ladder_score, dat$lw,
labels = dat$data$name,
xlab = paste("Happiness Score (", input$selected_year, ")", sep = ""),
ylab = "Spatially Lagged Happiness Score",
zero.policy = TRUE)
})
# LISA map
output$lisa_map <- renderTmap({
tmap_mode("view")
tm_shape(local_moran()$data) +
tm_fill(
col = "cluster_label",
palette = c(
"Insignificant" = "#ffffcc",
"Low-Low" = "blue",
"Low-High" = "#78c679",
"High-Low" = "#c2e699",
"High-High" = "red"
),
title = paste("LISA Cluster (", input$selected_year, ")", sep = ""),
style = "cat",
id = "admin",
popup.vars = c(
"Country" = "admin",
"Cluster Type" = "cluster_label",
"Happiness Score" = "ladder_score"
)
) +
tm_borders(alpha = 0.4) +
tm_layout(frame = FALSE, legend.outside = TRUE)
})
}
shinyApp(ui = ui, server = server)
```